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Abstract 



We investigate the phase diagram and the critical properties of the adsorbate 
. system sulphur/ruthenium(0001) in the coverage region | < < ^ by means 

I of Monte-Carlo simulations of a simple lattice gas model on a triangular lattice. 

■ The model contains only repulsive nearest and next-nearest neighbor interactions. 

On ■ Combining results obtained by using both Glauber and Kawasaki kinetics in the 

simulations we identify two tricritical points, three coexistence regions and a eutec- 
tic point in the phase diagram of the system. Our results agree with the findings of 
recent experimental work on S/Ru(0001). Furthermore we are able to add certain 
i details to the experimentally observed phase diagram. 

O 

CJ ■ PACS 0550 - Lattice theory and statistics; Ising problems 

PACS 6470 - Phase equilibria, phase transitions, and critical points 
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The comparison of experimental phase diagrams and of experimental data for the crit- 
ical properties of adsorbed layers with the results of computer simulations of appropriate 
lattice gas models can provide detailed information about the effective lateral interactions 
between the adsorbed particles. However, the reliability of this method depends on the 
validity of the lattice gas description for the system to be considered. Therefore, this 
point requires careful experimental tests. The inadequacy of the lattice gas description 
may be the reason for the quantitative disagreement between the simulated and the 
experimental p| phase diagrams of some systems. 

For the system S/Ru(0001), careful investigations [Q, ^, |[ have revealed that for 
coverages < 6 < | the S atoms are attached to well defined surface sites. For < < i 
only the the hep site is occupied while for higher coverages, | < 6 < |, one finds 
adsorbate atoms on the fee site as well. Thus, for coverages up to = | the description 
of the system by lattice gas models appears to be justified. Hence, for this system 
the comparison between Monte-Carlo (MC) simulations of such models and the rather 
complex experimentally observed phase diagram shown in Fig. |I| is well founded. Apart 
from the existence of five different commensurate structures the most interesting features 
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of the phase diagram of S/Ru(0001) are the striped domain wall phase and the coexistence 
region of two ordered structures which covers a large area in the interval | < 6 < |. To 
our knowledge, this is the first adsorbate system in which such a prominent coexistence 
region has been observed. In view of the complexity of this phase diagram it may seem 
questionable whether the experimental observations can be reproduced in detail by a 
reasonably simple lattice gas model. 

In this letter we shall demonstrate by means of MC simulations that, in fact, a surpris- 
ingly simple lattice gas model suffices for a quantitative description of the experimental 
phase diagram for coverages | < < | for which the coexistence region (region C) is 
observed. This area is bounded on the left by the homogeneous p(2 x 2) phase (B) and 
by the homogeneous (73 X 73) ph ase on the right (D) (Fig. 0). The ideal coverages 
(also referred to as critical coverages) of the two phases are © = | and © = | resp.. The 
coexistence region consists of homogeneous islands with p(2 x 2) and (v^ x ^/3) order. 

The lattice gas model with the smallest number of parameters which is able to produce 
these two structures is 

nn nnn 

As usual, Cj = 0, 1 denotes the occupation of the lattice site i, and the sums run over the 
N X N sites of the triangular lattice for nearest (nn) and next nearest neighbors (nnn). 

Positive coupling constants (fi for the first- and second-nearest neighbor interactions 
stabilize the p{2 x 2) structure at © = | (Refs. @, I^Dj while the existence of the (v^ x -\/3) 
structure requires that ^il^i < |, see Ref. The formation of p{2 x 2) islands at 
low coverages. Fig. 1, indicates that in a realistic lattice gas model for S/Ru(0001) a 
weak attractive third neighbor interaction has to be included. As we have checked this 
interaction is of little importance for coverages Q > ^. Therefore we ignore it in the 
present work. Similar models have been considered in previous numerical studies of 



adsorbate systems on square lattices and of systems on triangular lattices |T0|, [Tl|, |g, 
1^, For the parameter ratio (p2/fi = Jq^ ^ Monte-Carlo simulation of the model 



(|I]) has been performed by Glosli and Plischke 0, but a number of questions remained 
unresolved in their work. 

Commonly, simulations have been performed with the use of Glauber kinetics i.e. for 
constant chemical potential /x (Ref. fl^)- This requires that a term /i X] Q is added to the 
hamiltonian (0) in order to control the coverage © = -^J^iCi- Practically, with this tech- 
nique it turns out to be rather difficult to tune to a coverage that deviates significantly 
from the ideal coverages. In a Monte-Carlo run the system tends to jump between the 
competing ideal structures, and very long runs are necessary to achieve accurate statis- 
tics. Furthermore, this simulation technique does not directly mimic the experimental 
situation. The experiments are done at constant coverage and not at constant chemical 



potential. Thus it appears appropriate to use Kawasaki kinetics which conserves the 
coverage. 

To obtain a detailed picture of the phase diagram in the intermediate region between 
the ideal coverages © = |, | we combine results from both types of simulation techniques. 

As a measure for long range order of the p{2 x 2) type and of the (a/3 x a/3) type 
we use appropriate order parameters \l/p(2x2) and ^(^x^s) ^'^^ which explicit expressions 



can be found in Refs. [Q and [T^. In particular, these quantities allow the independent 
measurement of the fraction of these structures which constitute the coexistence region. 
To characterize the nature of the phase transitions, we calculate the susceptibilities 



Xp(2x2), X{v/3xV3) and the specific heat c, 



X = ^{('^n)-{^n)') (2) 
{{e%) - {sn)') . (3) 



Here e is the energy per site and ( • • • ) denotes the average over Monte-Carlo sweeps. 
The critical divergences of c and x in the infinite system appear as maxima for finite 
N. From the size dependence of these maxima the critical exponents a,7, z/ and the 
transition temperature Tc can be determined [l^, |l^ . 

In simulations at constant chemical potential the coverage fluctuates as a function 
of the temperature 

Ae:=^((0^)-(e^)^) . (4) 

The quantity A0 does not necessarily couple to the order parameter. However, it is 
clearly of interest for first order phase transitions: While the chemical potential and the 
temperature have to be equal in both phases the coverage jumps discontinuously at the 
transition. 

In our simulations we work with triangular lattices with periodic boundary conditions. 
The lattice sizes vary between 12 x 12 and 60 x 60 sites. Because of the different unit 
cells of the p{2 x 2) and and of the (v^ x ^/3) structures we have to change the linear 
dimension N in steps of 12 to avoid frustration by boundary effects. Each data point is 
obtained from 600 000 to 1 200 000 Monte-Carlo sweeps. 

The ratio of the order-disorder transition temperatures at the ideal coverages = |, | 
measured in the experiment is used to fix the coupling constant ip2 in units of ipi. As 
usual we set ipi = 1 and measure all energies in units of this nn coupling. This allows 
the direct comparison of results from both types of simulations. Agreement with the 
experiment is established with (f2 = 0.123. 

At the coverages = | both kinetics can be used for the determination of critical 
properties because the chemical potential is practically constant across the transition 
temperature. By comparison of the critical exponents with the results obtained from the 



experiment and with the values predicted by the Landau rules HlQl , |20| we are able to 
check our routines. Our finite size analysis reveals that Glauber kinetics and Kawasaki 
kinetics lead to identical results for both coverages: The maximum of the specific heat at 
= 1(1) has a size dependence c^'^^ ~ A^^ with f = 1.10 ± 0.1(0.50 ± 0.06), see Fig. |, 
while the susceptibility {xn"'^ ~ A^^) yields - = 1.95 ± 0.2(1.65 ± 0.2). These exponents 
are comparable with the exact values of the g = 4 (g = 3) Potts model - = 1(0.4), ^ = 
1.75(1.73) (Ref. [Q) and with the experimental values. Previous simulations of lattice 
gas models with p(2 x 2) and (y/S x 1/3) symmetries have yielded similar differences 
between the lattice gas exponents and the Potts exponents ll^, |ll|, ||, |12|, |13| . 

As we have discussed above, Kawasaki kinetics is the appropriate tool to analyze the 
ordered structures in the intermediate coverage region, | < < |, and the boundaries 
between these structures. For sufficiently low temperatures we find that particles added 
to the p{2 X 2) structure do not appear as point defects but as islands of (y/S x y/S) order. 
Due to the interfacial tension this minority phase forms hexagons in a sea of the ideal 
p{2 X 2) phase {majority phase). At ^ 0.292 ~ I (i + |) both domains occupy half 
the available area. Above = 0.292 the roles are interchanged, i.e. the denser (v^ x a/3) 



structure becomes the majority phase which surrounds islands of p{2 x 2) phase. The 
coexistence hne of the region III, Fig. (I), in which the p{2 x 2) and the (^3 x ^3) 
structure coexist, is determined by the disappearance of the order parameter measuring 
the minority domains. 

Clearly the boundaries of region III must end at the ideal coverages for T = 0. We 
are unable to confirm this for 6 < 0.27 and B > 0.31 because the area covered by the 
minority phase is too small. 

To avoid Fisher renormalization p2| we investigate the order-disorder transition of 
the majority phases, p{2 x 2) and (-\/3 x y^) resp., with Glauber kinetics. Paths of 
constant /i are shown in Fig. ^. As B approaches B = 0.292 the N-dependence of the 
finite-size peak in the specific heat c ~ N~ deviates significantly from the Potts values 
^ = 1, |, resp.. For yU = —0.50, label a, which corresponds to the ideal coverage B = 
the exponent is 1.10, see above. For ^ = —0.58 we find 1.33 (label b) and on the path 
fi = —0.68 we obtained 1.67 (label c). As B approaches B = 0.292 from higher coverages 
- increases from 0.5 for fi = —0.98, to 1.1 for fi = —0.74 (paths d - f). This behavior 
of ^ is indicative of the existence of tricritical points on the order to disorder transition 

lines. In the interval between these tricritical points Btr"* < B < the transition is a 
fluctuation induced first order transition. While the coverage remains close to the ideal 
coverages B = ;j, | resp. on lines of constant chemical potential which cross the phase 

boundary outside the interval Bt^'' < B < B^r"*, it changes rapidly on lines fi = const, 
that cross the phase boundary within this interval. On these last lines the coverage 
develops large fluctuations. This is seen in Figs. (^ and (^ where we plot AB versus 
kT. Evidently, for Btr'* < B < Btr"* the peaks in AB scale with a power of the system size. 
By contrast, AB is structureless in the vicinity of the ideal coverages, where the transition 
is of second order. The insets in Fig. ^ and Fig. ^ display the jump in coverage on lines 
of constant fi and give an estimate of the width of regions I, II in which the p{2 x 2) 
phase and the (v^ x 1/3) phase coexist with the lattice gas phase. From this point of 
view the maxima of the specific heat and of the susceptibilities obtained from simulations 
are to be considered as averages over these coexistence regions. The determination of 
the tricritical points from finite size results is difficult since the exponents ^, ^ increase 
smoothly when B is varied between the ideal coverages and B = 0.292. In general, the 
tricritical points are the lower (upper) end points of the coexistence regions I and II. 
We determine the boundaries of these coexistence regions by extrapolating the points of 
maximum curvature of the graphs B(T) shown in Fig. ^ and Fig. |^. This provides the 
estimate for the locations of the tricritical points in Fig. ^. 

The two upper boundaries of the coexistence regions I, II must meet the boundary of 
the coexistence region III in one point, i.e. at B = 0.292. Otherwise the region III would 
have a coexistence line with the lattice gas phase at which three phases coexist. This 
would contradict Gibbs' phase rule [jl9[ which, for a two component system with three 
phases, allows only for a triple point. This point is well known for binary alloys in 3d 
and is called eutectic point. We check the first order nature of the phase transition at 
the eutectic point with simulations for constant B. For symmetry reasons the chemical 
potential is constant on the line of constant coverage B = 0.292. Thus, in contrast 
to the general case B 7^ 0.292, calculations of critical exponents on this line are not 
affected by Fisher renormalization P^. At the eutectic point our simulations yield the 



following size dependence of the specific heat and the susceptibilities: c ~ A^^-^, Xp(2x2) ~ 
^^^^'^(Vsx Vs) ~ A^^'^^. For first order transitions these thermodynamic quantities are 
predicted to diverge as A^^ with the system size p3l. The reason for the discrepancy 



between this theoretical prediction and our numerical values for the exponents of the 
susceptibilities is visible in Fig. |. While the positions of the peaks of the specific heat, 
^(max)^ increase linearly with In(A^), the peaks of the susceptibilities, Xn'^^^ curve upwards 
as function of In(A^). This means that with our system sizes we have not yet reached the 
scaling regime. In fact, the data for Xp(2x2) and X{y^xV3) based on effective system 
sizes iVeff < Y because for 6 = 0.292 the p(2 x 2) and the {^/S x y/S) structures cover 
less than half the lattice. 

In summary, we have shown that the simple lattice gas model yields a surprisingly 
rich phase diagram in two dimensions which closely resembles phase diagrams commonly 
observed in three dimensional substitutional binary alloys. It accounts for all details of 
the experimental phase diagram in the coverage region ^ < O < |. This includes the in 
particular coexistence regions I and II. These regions and their structure could not be 
determined with certainty in the experiment although their existence had been inferred 
from the structure of the experimental phase diagram in the vicinity of = 0.292. 
Finally, we note that the topology of the phase diagram found in this study can be 
expected to occur more generally for systems with threefold symmetry provided the 
lateral interactions decrease with distance comparably fast as in our model. 
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Figure Captions 



Figure 1: Experimental phase diagram of the system S/Ru(0001) [Q. A: p{2 x 2) islands 
+ lattice gas. E': (y^ x ^/3) + hne defects. Regions A, E and E' are not considered in 
the simulation. 



Figure 2: Specific heat c for fixed coverage B = the curves correspond to the system 
sizes 24 x 24, 36 x 36, 48 x 48 and 60 x 60. The inset displays the linear fit through the 
maxima of c^. 



Figure 3: Calculated phase diagram. Solid lines: boundaries of different regions. Region I 
{II}: coexistence of the disordered phase with the p{2 x 2) {(-\/3 x v^)} phase. Region III: 
coexistence of p(2 x 2) and {VS x \/3). Tricritical points: Q[1^ = 0.256, fcTjt^^ = 0.173; 

(2) (2) 

e^;' = 0.324, kr;;' = 0.170 marked by •. Dotted lines: maxima of specific heat c for 
constant 9. Dash dotted lines: paths of = const.. 



Figure 4: Coverage fluctuations AG at first order phase transitions for three different 
values of /i corresponds to the labels (a - c) in Fig. @. Curves are shown for the same 
system sizes as in Fig. (^. 



Figure 5: Same as in Fig. ^ but for 6 > 0.292. 



Figure 6: Specific heat c (a); susceptibilities Xp(2x2), X(^x Vs) (b) the eutectic point. 
Insets: positions of the maxima of c, Xp(2x2) and X(y^xV3) functions of the lattice size 
A^. System sizes are in Fig. 
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